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ABSTRACT 

William Kyle Anderson, Doctor of Philosophy, 1986 

Major: Engineering, Department of Aerospace Engineering 

Title of Dissertation: Implicit Multigrid Algorithms for the Three- 

Dimensional Flux Split Euler Equations 

Directed by: David L. Whitfield 

Pages in Dissertation: 83. Words in Abstract: 232. 

ABSTRACT 

The full approximation scheme multigrid method is applied to 
several implicit flux-split algorithms for solving the three- 
dimensional Euler equations in a body fitted coordinate system. Each 
uses a variation cf approximate factorization and are implemented in a 
finite volume formulation. The algorithms are all vectorizable with 
little or no scalar computations required. The flux vectors are split 
into upwind components using both the splittings of Steger-Warming and ” 
Van Leer. Results comparing pressure distributions with experimental 
data using both splitting types are shown. The stability and smoothing 
rate of each of the schemes are examined using a Fourier analysis of the 
complete system of equations. Results are presented for three- 
dimensional subsonic, transonic, and supersonic flows which demonstrate 
substantially improved convergence rates with the multigrid algorithm. 
The influence of using both a V-cycle and a W-cycle on the convergence 
is examined. Using the multigrid method on both subsonic and transonic 
wing calculations, the final lift coefficient is obtained to within 0.1 
percent of its final value in as few as 15 cycles for a mesh with over 
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210,000 points* A spectral radius of 0.89 is achieved for both subsonic 
and transonic flow over the ONERA M6 wing while a spectral radius of 
0.83 is obtained for supersonic flow over an analytically defined 
forebody. Results compared with experiment for all cases show good 


agreement 



V 


TABLE OF CONTENTS 

ACKNOWLEDGMENTS ii 

ABSTRACT iii 

NOMENCLATURE vi 

LIST OF FIGURES 

LIST OF TABLES xii 

Chapter 

I. INTRODUCTION 1 

II. EULER SOLUTION METHOD 5 

2.1 Euler Equations in Generalized Coordinates ....5 

2.2 Flux Vector Splitting 7 

2.3 Baseline Solution Algorithm 15 

2.4 Boundary Conditions ••••••23 

2.5 Stability Analysis ..25 

III. MULTIGRID 30 

3. 1 General Algorithm ••••••••••• 30 

3.2 Algorithm for the Euler Equations 35 

IV. RESULTS 43 

4.1 ONERA M6 « • 43 

4.2 Analytic Forebody 61 

V. CONCLUSIONS 65 

APPENDICES 67 

A. Transformation to Generalized Coordinates ....67 

B. Splitting the Flux Vectors in 

Generalized Coordinates 71 

C. Restriction Operators 77 


BIBLIOGRAPHY 


81 



vi 


A,B,C 



A 


a 

CFL 

c 


F,G, H 


mass 


'energy 





if j 
J 




k i 


L 


NOMENCLATURE 

flux Jacobians , 9F/9Q, 9G/ 9Q, 9H/ 9 q 

A | A A . A A A 

flux Jacobians , 9F ± / 9Q, 9 G ± /dQ, Sh*/ 9Q 
matrix from similarity transformation 
speed of sound 

Courant-Friedrichs-Lewy number 

airfoil chord 

section lift coefficient 

pressure coefficient 

total energy per unit volume 

fluxes of mass, momentum, and energy 

mass component of flux vector 

energy component of flux vector 

grid level i 

identity matrix 

restriction operator used for transferring functions on 
grid i to grid i -1 

prolongation operator used for transferring functions on 
grid i -1 to grid i 

restriction operator for the residual 
/-I 

x and y components of a unit vector 
transformation Jacobian 

components of a unit vector normal to a 5 = constant face 
denotes £, n, or £ 
nonlinear operator 



vii 


M 

A 

N 

P 

P 




Q 


R 

A 

R 


r ,r ,r 
x y z 


S 


T 


T 


-1 


T 

t 


t 

x 


9 


t 

y 




z 


U,v,w 


U,V,W 


U 


0 


local Mach number in 5 direction 
implicit operator, also denotes Mach number 
unit vector normal to a k = constant face 
forcing function added to residual 
pressure 

approximate value of conserved variables on grid i 
conserved variables representing mass, momentum, and total 
energy per unit volume 

residual, also used as Riemann invariant 

A A 

unit vector normal to N and T 

A 

components of R in the x,y and z direction 
forcing function 
length of a cell face 

similarity transformation matrix whose columns are the 
right eigenvectors of the flux Jacobians, also used as a 
rotation matrix 

similarity transformation matrix whose rows are the left 
eigenvectors of the flux Jacobians 

A A 

unit vector normal to N and R 
time 

A 

components of T in the x,y, and z directions 
Cartesian velocities, v also denotes an eigenvector for 
Fourier analysis 

contravariant velocities, V also denotes the volume of a 
computational cell 
a constant vector 



v i 

x,y,z 

a 


a» 8, Y 


Y 

6 

X 

£/ h* ? 

A 

f • • • r Ag 

X 

u 

aq 

At 


Ax , Ay , Az 


3 

P 

I 


T 

T. 

1 

grad ( ft) 



correction on grid level i 
Cartesian coordinates 
angle of attack 

coefficients for stability analysis 
ratio of specific heats, Y= 1«4 
difference operator in x direction 
general curvilinear coordinates 
diagonal matrix of eigenvalues 
eigenvalues of Jacobian matrices 
amplification factor 
smoothing factor 

change in dependent variables, Q n+1 
time step 

change in x,y, and z 

denotes partial derivative 

density 

summation 

time 

relative truncation error 
gradient of ft 
magnitude of vector x 


Subscripts 

i,j,k cell indices 

x,y,z spatial derivative 


00 


denotes conditions at infinity 



N,N-1, ... 


discretization on grid 


Superscripts 

c denotes most current value of a quantity 

n denotes time level 

± denotes positive and negative flux and eigenvalue 

contributions; also denotes forward and backward spatial 
differencing or extrapolation 

( ) denotes quantities in locally orthogonal coordinates 

denotes fluxes and flux Jacobians in generalized 
coordinates; also denotes Fourier symbol 



X 


LIST OF FIGURES 

Figure Page 

1 • Variation of Steger -Warming split mass flux with 

Mach number • ...•••• . . 11 

2. Variation of Van Leer split mass flux with 

Mach number • 14 

3. Computational cell indexing 17 

4. Computational modules 

(a) 3-factor , spatially split scheme 19 

(b) 2-factor, eigenvalue split scheme*. 19 

(c) 2-factor, combination split scheme 19 

5. Vectorization of LU decomposition and back-substitution 

over multiple planes 20 

6. Stability analysis of three-dimensional approximate 
factorization schemes; M m = 0*8; a = 0°* 

(a) 3-factor, spatially split scheme 28 

(b) 2-factor, eigenvalue split scheme*.... 28 

(c) 2-factor, combination split scheme 28 

7. Multigrid V-cycle 39 

8 • Multigrid W-cycle 40 

9. Surface mesh for the ONERA M6 wing 

(a) 97 x 17 x 17 C-H mesh 44 

(b) 97 x 17 x 17 C-0 mesh 44 

10. Effect of multigrid on convergence; ONERA M6? 

193 x 33 x 33 C-H mesh; M = 0.84; a = 3.06°. 

CO 

(a) Residual history 45 


(b) Lift history 


45 



xi 

11. Comparison of convergence rate for the 3 schemes; 

ONERA M6 wing; 97 x 17 x 17 C-H mesh; 

M = 0.84; a = 3.06° 47 

00 

12. Comparison of convergence rate for a V-cycle and a 
W-cycle; ONERA M6 wing; 97 x 17 x 17 C-H mesh; 

= 0.84; a = 3.06° 49 

13. Upper surface variation of pressure coefficient; 

ONERA M6 wing; = 0.84; a = 3.06°. 

(a) 193 x 33 x 33 C-0 mesh 52 

(b) 193 x 33 x 33 C-H mesh 52 

14. Pressure coefficient contours on upper surface; 

ONERA M6 wing: = 0.84; a = 3.06° 53 

15. Comparison of experiment and inviscid calculations 
using Van Leer splittings; ONERA M6 wing; 

M = 0.84; a = 3.06° 54 

00 

16. Comparison of experiment and inviscid calculations 
using Steger-Warming splittings; ONERA M6 wing; 

M = 0.84; a = 3.06° 55 

00 

17. Multigrid convergence for C-H and C-0 meshes; ONERA 

M6 wing; = 0.699; a = 3.06°, V-cycle .........57 

18. Comparison of experiment and inviscid calculations 
using Van Leer splittings; ONERA M6 wing; 

M = 0.699; a = 3.06° 59 

CO 

19. Comparison of experiment and inviscid calculations 
using Steger-Warming splittings, ONERA M6 wing; 

M = 0.699; a = 3.06° 60 

00 

20. Grid for analytic forebody 62 

21 • Convergence history for analytic forebody; 

M =1.7; a = 0° 63 

CO 

22. Calculated pressure coefficient comparison with experiment 

at forebody symmetry plane; = 1.7; a = 0° 64 

23. Control volumes for restriction of dependent 

variables and residual..... 78 



XXI 


LIST OF TABLES 

Table Page 

1 Operation counts for solving the implicit operators 24 

2 Summary of results for ONERA M6 wing 

M = 0.84; a = 3.06° 50 

00 

3 Summary of results for ONERA M6 wing 

M = 0.699; a = 3.06° 

00 


58 



1 


Chapter I 
INTRODUCTION 

Upwind difference schemes for solving the Euler equations are 
becoming increasingly popular for several reasons. The time-dependent 
Euler equations form a system of hyperbolic equations and upwind 
differencing models the characteristic nature of the equations in that 
information at each grid point is obtained from directions dictated by 
characteristic theory. Some of the methods include the A-method, 1 the 
split-coefficient method, ^ flux-vector splitting, ^#4, 5 and f in- 
difference splitting. 6 These methods can be classified, in general, as 
upwind methods and have the advantage of being naturally dissipative. 
Separate spatial dissipation terms, such as are generally required in a 
central difference method to overcome oscillations or instabilities 
arising in regions of strongly varying gradients, need not be added. 

While the A-method and the split-coefficient method closely mimic 
the method of characteristics, they are both applied to the 
nonconservative form of the equations and consequently require the use 
of shock-fitting techniques to obtain the correct location and strength 
of shocks in transonic flows. Use of the conservation-law form allows 
shock waves to be captured as weak solutions to the governing equations 
and circumvents the difficulty in applying shock-fitting techniques to 
arbitrary flows. Both the flux-difference-splitting and flux-vector- 
splitting methods can be applied to the conservation-law form. 

The particular upwind method used in the current work is the flux 
vector splitting method in which the flux vectors are split into for- 
ward and backward contributions based on an eigenvalue decomposition and 
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differenced accordingly. The splittings investigated include those of 
Steger-Warming 3 ' 4 ' 7 and Van Leer. 5 ' 8 ' 9 The advantages of flux splitting 
are obtained at the cost of increased computational work in comparison 
to unsplit methods, since two sets of fluxes are computed for each 
coordinate direction and implicit schemes require two sets of flux 
Jacobians (e.g. 8 f + / 3Q and 8 f"/ 3Q) for consistent linearization of the 
fluxes. In addition, the split fluxes and flux Jacobians are also 
generally more complicated than the unsplit terms owing to the branching 
involved with eigenvalue sign changes. 

In order to offset the additional computational work of the upwind 
methods, it is highly desirable to accelerate the convergence rate, 
especially when only steady-state solutions are sought; the objective is 
to reduce the computer time required while still maintaining the high 
level of robustness and accuracy attained from upwind differencing. 
Accelerating the convergence rate becomes increasingly important as the 
mesh is refined since the log of the spectral radius for single grid 
methods generally varies linearly with the mesh size, making 
computations on very fine meshes impractical. 

One method which has been successful in accelerating the 
convergence rate of elliptic problems, attaining a spectral radius 
independent of the mesh spacing, is the multigrid method. 10 ' 11 Although 
most of the existing theory on multigrid methods pertains specifically 

in n n 

to elliptic equations, it has been shown in several references that 

the multigrid method can greatly accelerate the convergence rate of 
numerical schemes used for solving the Euler equations. 
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One of the earliest applications of multiple grids in solving the 
Euler equations, was presented by Ni who used coarse grids to rapidly 
propagate corrections throughout the domain. 12 His original idea was 
first incorporated into a one-step Lax-Wendroff method and later 

1 3 

extended for use into predictor-corrector type methods by Johnson. 
Johnson 1 ^' 18 and chi ma and Johnson 1 ^ subsequently used the method to 
calculate both inviscid and viscous flows over several two dimensional 
geometries. In 1984, Mulder applied a linear multigrid scheme to the 
Euler equations in two space dimensions using upwind differencing to 
calculate flow over a circular arc and for a weakly barred galaxy. 1 ^ 
Jesperson also used upwind differencing in two spatial dimensions to 
calculate flow over airfoils. 18 In this approach, the Euler equations 
were solved by Newton iteration where the linear system arising at each 
step was solved using multigrid. One of the first uses of the nonlinear 
multigrid method in accelerating the convergence rate for both the two 
and three dimensional Euler equations was reported by Jameson who used 
central differencing in a four stage Runge-Kutta algorithm to advance 
the solution. 1 ^ 9 20 In two dimensions, recent work by Jameson and Yoon 
also used central differencing and incorporated the multigrid algorithm 
into some implicit schemes with good success. 21 ' 22 

The purpose of the current investigation is to combine the full 
approximation scheme (FAS) multigrid method with flux vector splitting 
to obtain efficient solutions to the Euler equations in three- 
dimensions. The full approximation scheme for a general nonlinear 
problem is discussed as well as its implementation for the Euler 
equations. Both V and W type cycling strategies are also 
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investigated. Several smoothing algorithms are given involving mostly 
vectorizable computations on the VPS-32 supercomputer at NASA Langley. 
In addition, both the splittings of Steger-Warming and Van Leer are 
considered for splitting the flux vectors into upwind components. 
Numerical pressure distributions are compared with available experiment 
for subsonic and transonic flow over the ONERA M6 wing and supersonic 
flow over an analytically defined forebody. 
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Chapter II 

EULER SOLUTION METHOD 
2.1 Euler Equations in Generalized Coordinates 

The governing equations are the time-dependent equations of ideal 
gas dynamics/ i.e., the Euler equations, which express the conservation 
of mass, momentum, and energy for an inviscid nonconducting gas in the 
absence of external forces* The conservation form of the equations in 
generalized coordinates is given by 


8Q 8f 8G 8h 
3t + 3£ + 3n + 3c " 


0 


where 



F 


J 


PU 

pUu + 5 p 

X 

pUv + 5 P 

y 

pUw + £ p 
z 

(e + p)U 


( 2 . 1 ) 


(2.2) 


(2.3) 


I pV ' 

pvu + n p 
x 

pvv + n y p , 

pvw + n p 

z 

(e + p)V 


(2.4) 
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H 


_ 1 _ 

J 


PW 

pWu + £^p 
pWv + £ p 

y 

pWw + X> P 
z 

(e + p)W 


(2.5) 


The pressure p is related to the conserved variables through the ideal 
gas law 


P = (Y - D [e - p (u 2 + v 2 + w 2 )/2] 


( 2 . 6 ) 


The equations have been generalized from Cartesian coordinates using a 
steady transformation of the type 


£ = £(x, y, z), n= n(x, y, z) , c = C(x f y, z) , x = t 


(2.7) 


where the contravariant velocity components are 

U = 5u + 5'v + 5w (2.8a) 

x y z 


v = n u + n v + n w 
x y z 


(2.8b) 


W = C x u + ^v + ^w (2.8c) 

The transformation to generalized coordinates is given in Appendix A. 

The equations, while written in generalized coordinates, are used 
in a finite-volume formulation. Equation (2.1) can be interpreted as 
describing the balance of mass, momentum, and energy over an arbitrary 
control volume. In this connection, the vectors grad(£)/Jf grad( r])/J, 
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and grad(£)/J represent directed areas of cell interfaces in the 
5, ri/ and £ directions and the Jacobian J represents the inverse of the 
cell volume. Likewise, the quantities pU/J, pV/J, and pW/J represent 
the mass flux crossing the cell interfaces in the £, rjr and ? 
directions • 

2.2. Flux Vector Splitting 

The upwind differencing in the present work is effected through 
the technique of flux-vector splitting. The generalized 

A A A 

fluxes F, G, and H are split into forward and backward contributions 
according to the signs of the eigenvalues of the Jacobian matrices and 
differenced accordingly. For example, the flux in the £ direction can 
be differenced as 

6 F - F + + 6* f“ (2.9) 

*9 ^ S 

A + , A - 
since F has all non-negative eigenvalues and F has all non-positive 

eigenvalues. For the current study, two methods of splitting the flux 

vectors into upwind components are considered. Although the details 

of each method can be found in references 3, 4, 5, 7, and 9 both 

methods are briefly discussed below. 

The first method presented is the technique outlined by Steger 

and Warming in reference 3. Since the flux vectors are homogeneous 

A 

functions of degree one in Q they can be expressed in terms of their 
Jacobian matrices. Considering the flux vector in the ^-direction for 

A 

example, F can be written as 
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* A A Cj-ri * 

F = AQ = — ^ Q 

8Q 


( 2 . 10 ) 


Using a similarity transformation, equation (2.10) can be rewritten as 


A A A 


-1 


F = AQ = TAT Q (2.11 ) 

The matrix A is a diagonal matrix composed of the eigenvalues 

A 

of A and is given by 


A = 


0 

X 2 

0 

0 

0 


0 

0 

x 3 

0 

0 


0 

0 

0 


0 

0 

0 

0 

X. 


( 2 . 12 ) 


X. _ _ =U = £u + 5v+5w 
1,2,3 x y 


where X^ = U + |grad(5)|a 
X5 = U - | grad ( £) |a 


(2.13) 


The eigenvalues can then be decomposed into non-negative and non- 
positive components 


x. = x + + x. 

11 1 

where 


X. 

i 


± 



(2.14) 


(2.15) 
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Similarly, the eigenvalue matrix A can be decomposed into 
A = A + + A" (2.16) 


where A + is made up of the non-negative contributions \^ + and A is 
constructed of the non-positive contributions \ . This splitting of 
the eigenvalue matrix, combined with equation (2.11) allows the flux 

A 

vector F to be rewritten as 

F = T(A + + A*)T _1 Q = (A + + A")Q = F + + F" (2.17) 


The flux vector F has three distinct eigenvalues given by (2.13) 
and can therefore be written as a sum of three subvectors, each of which 

7 

has a distinct eigenvalue as a coefficient. 


where 


F = F + F + F 
12 3 


(2.18) 


= A. 


P 

pu 

Y-1 ] Pv 

i JY I pw 

f P , 2 2 2 X 

| — (u + V + w ) 


(2.19) 


2,3 


= A 


PU ± pa E 

A* 

pv ± pa E 

aY 


4,5 J2y I pw ± pa 5 


e + p ± 


paU 


grad ( £) 


(2.20) 
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and the direction cosines of the directed interface in the ^-direction 
are 

i = /k ad <s> 

\ - 5 y /|grad<5) 
i = 5 /I grad ( 5) 

Z Z 1 

The forward and backward flux vectors F and F are formed from 

equations (2.18), (2.19), and (2.20) by inserting 

X. = X. + and X. = X. ", respectively. It should be noted that for 

li ii 

supersonic and sonic flow, in the ^-direction, i.e., |m^| = |u/a| > 1 r 

where u = U/|grad(5)| represents the velocity normal to a 5 = constant 
face, the fluxes in this direction become 


(2.21a) 

(2.21b) 

(2.21c) 



( 2 . 22 ) 


The split fluxes in the other two directions are easily obtained by 
interchanging t] or £ in place of £. 

The fluxes split in this manner above are not continuously 
differentiable at zeros of the eigenvalues (i.e. sonic and stagnation 
points) .23 This is illustrated in figure 1, where the split mass flux 
contributions for the one-dimensional Euler equations, non- 
dimensionalized by division through pa, are shown as a function of the 
Mach number. The gradient discontinuities in the split fluxes are 
evident as the eigenvalues pass through zero. The lack of 


differentiability of the split fluxes has been shown in some cases to 
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cause small oscillations or glitches at sonic points but which are 
rarely noticeable for most aerodynamic applications. 

A, A 

It should also be noted that the Jacobian matrices of F + and F~ 

which are required for proper linearization for an implicit scheme do 

not have the same eigenvalues as A + and A~ defined in equation 

(2. 17). 23 However, the Jacobian matrices of F and F do have the same 

sign as A + and A so that upwind differencing the spatial derivatives 

remains appropriate.^ Although they are easier to form, the use 

A + 

of A in implicit schemes instead of the correct linearizations A has 

O nA 

been shown in many cases to cause severe time step limitations. ' 

In 1982, a new method of splitting the flux vector was proposed by 
Van Leer. 5 The approach taken in the derivation was to split the fluxes 
so that the forward and backward flux contributions transitioned 
smoothly at eigenvalue sign changes, i.e., near sonic and stagnation 

points. Just as for the Steger-Warming splitting, it was required that 

3f^ . . 3F 

the Jacobian matrices — — have non-negative eigenvalues and — — have 

3Q 3Q 

non-positive eigenvalues so that upwind differencing could be used for 

the spatial derivatives. In addition it was required that both 

Jacobians have one zero eigenvalue for subsonic Mach numbers which leads 

to steady transonic shock structures with only two transition zones. 5 

In practice, when second-order spatial differencing is used, shocks with 

only one interior zone are usually obtained. 9 This feature is not 

observed with the Steger-Warming flux splitting. 

The three-dimensional splittings of Van Leer were originally given 

for Cartesian coordinates. The extension to generalized coordinates is 

given in Appendix B with the resulting split fluxes given below. Only 
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the splitting for the flux in the 5-direction is given, as the others 

A 

can be obtained similarly. The flux vector F is split according to the 
contra variant Mach number in the ^-direction, defined above 
as = u/a. For supersonic flow, 


+ 

F = F, 


F = 0, 

A A A 

F~ = F, F = 0, 

and for subsonic flow, | M | < 1 


> + 1 
M < - 1 


r>± _ grad (5) 
F — j — 


r f ± 

mass 

"mass 1 '* 1 - 5 * * «1 


r± 

'me 
'me 

B* 

energy 


f x [5 (- u ± 2a)/ y + v] 
mass y ' 


f mass U z ( " U ± 2a) / Y + w] 


where 

f ± 

mass 


= ± pa(M^t 1) / 4 


(2.23) 


(2.24a) 


(2.24b) 


f* = f f aee l{-(Y - 1)u 2 ± 2(y - 1 )ua 

energy mass 


+ 2a 2 }/(y 2 - 1) + (u 2 + v 2 + w 2 )/2] 


(2.24c) 


For forming F*, are given by equation (2.21) and u is the 

velocity normal to a 5 = constant face. The fluxes in the other two 
directions are easily formed by interchanging £ with either n or 
In figure 2, the non-dimensionalized mass flux, using the Van Leer 
splitting, is shown as a function of Mach number for the one-dimension- 
al Euler equations. The split fluxes are continuously differentiable at 
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sonic and stagnation points; the improvement over the Steger-Warming 
splitting is apparent. 

2.3 Baseline Solution Algorithm 

The baseline algorithm for updating the steady Euler equations 
stems from a backward Euler time integration of the unsteady equations 

Q 

which yields 


[I + At( 6^A + +6^A") + At ( 6~B + +<S*B~) + At ( <5~C + + 6*C - ) ] AQ = -AtR n 


where the residual at time level n is given by 

R n = 6“ F + + 6* F" + 5" G + + 6 + G - + 6" H + + 6 + H 

5 S n T1 C C 


(2.25) 


The split-flux differences in equation (2.25) are implemented as a 

flux balance across a cell, corresponding to MUSCL type differencing 

25 

(Monotone Upstream-centered Schemes for Conservation Laws). For 
example, the flux balance in the ^-direction across a cell centered at 
point i,j,k can be written as 

6“ F + + 6+ F“ = [F + (Q~) + F-(Q + )]. +1/2 

- (F + (Q“) + F "(Q + )l i _ 1/2 (2.26) 

- 

The notation F (Q denotes the forward flux evaluated using the 

metric terms at the cell interface i+1/2 and the conserved state 


variables on the upwind side of the interface, obtained by a fully- 
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upwind second-order state variable interpolation: 


e "i+1/2 - u5 <%. - °' 5 
e i+1/2 " 1,5 °i+1 ~ °* 5 ®i+2 


As seen in figure 3, Q? . , denotes the average value of Q in the cell 

1 r 3 t K 

centered on at time t n ; for simplicity, wherever the script 

notation is i,j,k, or n, it is most often dropped. 

In equation (2.25), if second order differencing is used on both 
sides of the equation, Newton iteration for the steady Euler equations 
is obtained as At tends to infinity. The solution however requires the 
solution of a large banded block matrix at each step which is generally 
not feasible due to the amount of operations required to invert the 
system. Even if the differencing on the left hand side of the equation 
is reduced to first order, which would not effect the second order 
accuracy of the final solution, the resulting system of equations 
usually remains uneconomical to solve. Therefore, the solution is 
obtained using approximate factorization, which splits the implicit 
operator into a sequence of easily invertible equations. 

When using flux vector splitting, there are numerous ways of 

3 

factoring the implicit operator into a sequence of simpler operators. 
For the results shown below, three ways of factoring are considered. 
Each of the schemes uses first order spatial differencing on the 
implicit side of the equation while second order differencing is 
maintained for the residual calculations. Since the steady state does 
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not depend on the differencing of the implicit side, the final steady- 

state result will be spatially second order accurate. All the schemes 

employ simple explicit boundary conditions* Since only steady-state 

* 

solutions are sought, each cell is advanced at its own time step 
corresponding to a given CFL number defined by 

CFL = At[ |u| + |v| + |w|+a( |grad( 5) | + |grad( n) | + |grad( C) | ) ] (2.28) 

The first scheme considered is a spatially-split algorithm given by 


[I+At( 6~A*+6tA ) ] [I+At( S“b + +6 + b“)] 

5 Z n n 

[I+At( 6~C + +6*C~) ] AQ =-AtR n 


(2.29) 


The computational module for the implicit side of equation (2.29) is 
shown in figure 4a for the £ sweep. Since the solution at each point is 
directly coupled to the two neighboring points, the scheme requires the 
solution of a system of block tridiagonals. Similarly, the other two 
factors also require a block tridiagonal inversion. This scheme has the 
advantage however of being completely vectorizable and viscous effects 
can be easily incorporated into the left hand side. Since the speed of 
the VPS-32 is much faster for long vector lengths than it is for short 
ones, the computations in the present implementation take advantage of 
the large memory available on the VPS-32 computer and solve the block 
matrix equations over multiple planes simultaneously, yielding longer 
vector lengths and faster processing rates. As seen in figure 5, the 



(a) 3-factor spatially-split 


O Time level n 

• Time level n + 1 

© Time level n + 1 

(previously calculated) 





Positive eigenvalue 
sweep 


(b) 2-factor eigenvalue split 



^ E, /C + sweep 


Figure 4. Computational modules 
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block tridiagonal matrices can be solved with vector lengths 
corresponding to the number of lines in a plane times the number of 
planes taken. The residual calculations, on the other hand, can be made 
with vector lengths corresponding to the number of points in the grid. 

To decompose the implicit operator into lower and upper matrices (LU 
decomposition) and perform the back substitutions requires 695 
multiplications and additions per factor resulting in a total of 2085 
operations for each sweep through the grid. 

The second method considered to factor the left hand side of 
equation (2.25) is a two-factor method in which the implicit operator is 
split such that one operator contains the Jacobians with all positive 
eigenvalues and the other operator contains the Jacobians with all 
negative eigenvalues. The scheme can be written as 


[I+At( 6 "a + +6"b + +6 C + )] 

5 n C 


[I+At( 6^A - +5*B~+6*C~) ] AQ = -AtR n 


(2.30) 


This scheme only requires the solution of block lower triangular 
equations. Each factor is solved by starting at one corner of the grid, 
solving for each point by marching across the field to the opposite 
corner. From the computational module shown in figure 4b for the first 
factor in equation (2.30), the solution at the center node requires that 
the solution of each of the points behind it be previously obtained. 
These terms are taken to the right hand side of the equation and added 
to the residual. A similar procedure is carried out for the second 
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factor. For the present implementation of this scheme, only 270 
operations are required for each factor to invert the left hand side and 
perform the back substitutions and most of the algorithm is 
vectorizable. This scheme requires that a 5 x 5 matrix be inverted at 
each point in the grid which can be carried out with vector lengths 
corresponding to the number of points in the grid. However, since the 
solution of each plane requires that the solution of the previous plane 
be known, the back substitution cannot be performed over multiple 
planes. The maximum vector lengths in the back substitution process 
corresponds to the number of points in a plane. The only scalar 
computations correspond to back substitution along a line requiring 
roughly a third of the total operations. Although the scheme requires 
only about 25 percent of the operations required for the spatially split 
scheme, the scalar computations degrade the processing rate 
significantly so that the overall processing rate on the VPS-32 is about 
twice as slow as for the spatially split scheme. 

The last scheme considered is another two factor scheme which is 
spatially split in two directions, with the third direction split 
according to the sign of its eigenvalues. The resulting scheme, which 
is referred to as combination splitting, is given by 

[I+At( 6> + +6tA _ +6~C + )] [I+At( 5-B + +6 + B"+6tc“)] AQ = -AtR n (2.31) 

£f C C n n 5 

and the computational module for the first factor is shown in figure 
4c. This scheme also requires the solution of block tridiagonal 
systems, requiring about 695 multiplications and divisions for each 
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factor. This scheme is completely vectorizable; however, as in the 
previous two-factor scheme, the solution of each plane requires that the 
solution of the previous plane be known thereby eliminating the 
possibility of extending the vector operations over several planes. The 
result is that even though this scheme requires only two-thirds of the 
operations of the 3-factor scheme, the computational rate is actually 
degraded by about ten percent. A summary of operations required to 
solve the left hand side for each of the three schemes is given in table 
1 . 


2.4 Boundary Conditions 

The boundary conditions for the solutions presented below are 
applied explicitly. On the body, the normal velocity is set to zero 
while the pressure and density are determined by extrapolation from the 
interior. In the farfield, for subsonic flow, the velocity normal to 
the boundary and the speed of sound are obtained from two locally one- 
dimensional Riemann invariants given by 

R* = u ± (2.32) 

These are considered constant along characteristics defined by 

(^•) ± = u ± a (2.33) 

For subsonic conditions at the boundary, R~ can be evaluated locally 
from conditions outside the computational domain, and R + locally from 
the interior of the domain. The two Riemann invariants can be added and 


subtracted to determine a local normal velocity and speed of sound at 



Scheme 

Operations required per factor and vector length 

LU decomposition 

Additional right-hand terms* 

Backsubstltutlon 

3-factor spatially split 

550 (multiplane) 

• 0 

1A5 (multiplane) 

2-factor combination split 

550 (multlplane) + 

50 (plane) 

1A5 (single plane) 


Scheme 

Operations required per factor and vector length 

Additional right-hand terms* 

LU decomposition** 

(LHS)-l-(RHS) 

2-factor eigenvalue split 

50 (plane) 
50 (line) 

50 (scalar) 

75 (multiplane)"** 

A5 (scalar) 


"Terms required on right-hand side ( RIIS ) In addition to residual 
**LU decomposition of a 5 x 5 matrix only 

*length corresponds to the number of lines In a plane times the number of planes 
++ 1 ength corresponds to the number of points In a plane times the number of planes 


Table 1. Operation counts for solving the implicit operators 
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the boundary respectively. Depending on the sign of the normal velocity 
(inflow or outflow where inflow at the boundary corresponds 
to u < 0)/ the entropy and tangential velocities extrapolated from the 
exterior or interior of the domain are used with the speed of sound to 
determine the density and pressure on the boundary. The Cartesian 
velocities are obtained by decomposing the velocity tangential and 
normal to the boundary. 

For supersonic free-stream conditions along inflow boundaries, 
quantities are extrapolated from the exterior; along outflow boundaries, 
quantities are extrapolated from the interior of the computational 
domain. 

2.5 Stability Analysis 

In order to examine the stability characteristics of the three- 

dimensional approximate-factorization algorithms considered above, a 

Fourier analysis is conducted on the complete system of equations in 
• ft 

Cartesian coordinates. Because of the mixed signs of the 

eigenvalues of the Euler equations and the fact that the three- 
dimensional Euler equations cannot be diagonalized to yield a system of 
convection equations, stability analysis of the scalar convection 
equation is not sufficient to determine stability properties of the 
three schemes. The complete system of equations can be written as 

MAQ = -L =-AtR n (2.34) 


where for Cartesian coordinates 
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R n = 6 _ F + + 6 + F _ + <5 G + + 6 + G~+ 6~H + + 6 + H (2.35) 

x x y y z z 


and M is an implicit operator corresponding to the scheme considered. 

Linearizing the residual R n as 

_ n + r” ^n — ~4* _n _+ ~n 

R = A0Q + A0Q + B6Q 
x x y 

(2.36) 


+ B~6 + Q n + C + 6~Q n + C 6V 
y z z 

and assuming that the Jacobians are locally constant, the stability can 
be analyzed by letting 


^n .n TT iBx i7V i oe 
Q = A U e e ,J e 
o 


(2.37) 


where U Q is an initial constant vector. Upon substitution into equation 
(2.34) and dividing out the common factors, the generalized eigenvalue 
problem for A, which is the vector of amplification factors, can be 
obtained 


A A A 

(M-L)v = MAv (2.38) 

A A 

where M and L are the Fourier symbols of M and L, respectively. The 
stability characteristics are determined by cycling through a fixed 
number of each of the spatial frequencies, in this case sixteen 
frequencies, in the range 


0 < 3Ax, yAy, ctAz < 2 tt 
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for a series of CFL numbers between 0.1 and 50. The generalized 
eigenvalue problem is solved each time using a routine from the 
International Mathematics and Statistics Library (IMSL). 2 ^ Each time, 
the maximum eigenvalue, average eigenvalue, and the smoothing factor are 
determined, where the smoothing factor, defined as 

_ max {| A| } 

tt/ 2 < max(3Ax, yby, cxAz) < 3 tt/2 (2.39) 

corresponds to the damping of the high frequencies and serves as an 

indication of how effectively the multigrid procedure can accelerate 

convergence for a given scheme. 

Results are shown using the Van Leer splitting for each of the 
schemes given above. Identical cases were run using the Steger-Warming 
splittings with little change in the results. Each result was obtained 
by using first order differencing on the implicit side of the equation 
and fully-upwind, second-order differencing for the residual 
computations. All the calculations assume Cartesian coordinates, a Mach 
number of 0.8, and zero degrees yaw and angle of attack. 

The average eigenvalue, the smoothing factor, and the maximum 
eigenvalue is shown in figure 6 for the three schemes given above. For 
the 3-factor, spatially split scheme, shown in figure 6a, the maximum 
eigenvalue indicates that this scheme is conditionally stable with a 
maximum CFL number of approximately 20. The minimum smoothing factor 
occurs at a CFL of about 5 which is somewhat less than the CFL number 
where the maximum eigenvalue is lowest. In contrast to the spatially 
split algorithm, both of the 2-factor schemes shown in figure 6b and 6c 
appear to be stable for all CFL numbers considered. The maximum 
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eigenvalues and smoothing factors also exhibit less sensitivity to the 
CFL number than the 3-factor scheme with minimum smoothing factors also 
occurring at a CFL number of about 5. However, the 2-factor eigenvalue 
split scheme has a somewhat higher smoothing factor than the other two 
schemes which are therefore more appropriate for multigrid applications. 
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Chapter III 
MULTIGRID 

3*1 General Algorithm 

The multigrid method used in the current study is the full 

1 0 28 29 30 

approximation scheme (FAS) which appears in many references ''' 
and is summarized below. It is most easily understood by first 
considering the solution of a general nonlinear equation 

L(Q) = S (3,1) 

Equation (3.1) is to be solved numerically by dividing the domain 
into discrete cells yielding a system of equations to be solved 
simultaneously at each point as 


w - s » <3 - 2) 

where is the exact solution to the discretized system and L^ is the 
discrete analog of the operator L. If initial conditions are close 
enough to the final solution, equation (3.2) could be solved iteratively 
using Newton iteration. This approach however may be prohibitively 
expensive if the number of unknowns is large as typically occurs in 
multi-dimensional problems. Many other iterative schemes have therefore 
been devised which require significantly fewer operations. After a few 
iterations however, these methods generally exhibit a slow convergence 
rate, reducing the residuals by a very small amount each time.^ 9 


The 
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reason for the slow asymptotic convergence rate is inadequate damping of 
the low frequency errors • 1 1 

The multigrid method efficiently damps the low frequency errors 
using a sequence of grids Gq, , . . . G N where G N denotes the finest 
grid from which successively coarser grids can be formed by deleting 
every other mesh line. In this context, the high frequency error 
components on a given grid are those which cannot be resolved on the 
next coarser mesh due to the increased grid spacing. If an iterative 
method is chosen which quickly damps the high frequency errors on a 
given grid, the remaining errors will be the lower frequency smooth 
components. The sequence of coarser grids can then be used in 
accelerating the convergence rate on the finest grid by reducing the 
remaining low frequency errors since some of these same frequencies 
appear as high frequency errors on a coarser grid. Therefore, the 
errors on the fine grid which are responsible for slow convergence are 
quickly damped using the coarser grids where the computations are 
relatively cheap. 

In order to use the coarser grids, it is necessary to obtain an 
equation on the fine mesh which can be accurately represented by the 
coarser mesh. It is important to first realize that neither the high 
frequency solution or error components on the fine grid can generally be 
resolved on a coarser grid. The high frequency errors however can be 
sufficiently damped on a fine grid using a variety of iterative schemes 
so that the remaining errors will' mainly be composed of smooth low 
frequency components which can be adequately represented on coarser 
meshes. For this reason, it is necessary to obtain an equation on the 
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fine mesh in terms of the errors* 

When solving in an iterative fashion, equation (3.2) is solved 
approximately at each step as 


l n(1 C n) - % + % 


(3.3) 


where q c N is the most current approximation to Qjg and R N is the residual 
which will only be zero when q^ = and hence the exact solution is 
obtained. Subtracting (3.3) from (3.2) yields an equation on the finest 
grid in terms of the residual 


W - - -"n 


If the high frequency errors have been previously smoothed, then the 
fine grid residual equation (3.4) can be adequately approximated on a 
coarser mesh by 


N-1 


<W = i T 1 ( A ) + Vi CV> 


(3.5) 


N-1 ^N-1 

Where t and I are restriction operators for transferring the 
N N 

dependent variables and the residual from the fine grid to the coarse 
N-1 c . . 

grid. Here, 1^ q^ serves as an initial approximation to the solution 
on the coarse mesh whereas is the exact solution which is the sum 

of the initial approximation and a correction. Since the full 
solution is computed and stored on each grid level as opposed to only 
the corrections, this is referred to as the full approximation scheme 


(FAS) 
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On the coarser grid, equation (3.5) could possibly be solved 
exactly using a variety of numerical techniques to obtain Q N-1 from 
which the correction can be formed as 


V N-1 “ Qn-1 


- I, 


N-1 c 
N 


. (3.6) 


This can then be passed up to the fine grid and used as a correction to 
q N c which is replaced by its previous value plus the correction 




+ I 


N 

'N-1 


N-1 


(3.7) 


This process yields a simple FAS two level algorithm where the 
operations on the coarse grid (equations 3. 5-3.7) which are used to 
update the fine grid solution are termed the coarse grid correction. 
Often however, the exact solution of (3.5) can still be quite expensive 
to obtain. Also, since the correction on the coarse grid serves only as 
an approximation to the fine grid correction, the exact solution of 
(3.5) is not required. Therefore, instead of solving (3.5) to 
completion, several iterations can be carried out to get a reasonable 
approximation to . After each iteration of equation (3.5), the 

c 

equation satisfied by is given by 

Vi <•£-, > - 5" 1 <-v + Vi * Vi (3 - 8> 


which differs from the solution of (3.5) only by the residual term which 

Q 

will be zero when 0^=0 .If the errors are smooth, then 

N-1 TJ-1 
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subtraction of (3*8) from (3.5) yields an equation which can be well 
represented on yet a coarser mesh, Writing this equation on G N-2 

yields 


L (Q ) = L (I N_2 q C ) 
N-2 V N-2 N-2 N-1 4 N-1 


+ 3:? <-Vi> 


(3.9) 


where equation (3.8) defines • The solution to (3.9) may be solved 

for exactly, approximated by several iterations, or alternately, by 
introducing more coarse grid levels. On all coarse grids, one or more 
FAS cycles (smoothing followed by coarse grid correction) are done. In 
this manner, each of the coarse meshes is used to obtain a correction 
for the solution on the next finer mesh. Since only the equations for 
smooth error components may be represented well on coarser grids, it is 
important to only pass corrections from a coarse grid up to the next 

o Q 

finer one and not the full solution. As a result of using the coarser 
meshes, a fast convergence rate may be obtained since the low frequency 
error components on the fine grid are quickly damped on the coarser 
meshes. 

Note that equation (3.5) can be recast using equation (3.3) as 


l n-i (2 n-i 5 


N-1 


+ T . 


N-1 


= P 


N-1 


(3.10) 


where 


S = I* 5-1 s 
N-1 N N 


_ , N-1 C. "N-1 /T , C., 

T N-1~ L N-1 (I N q N “ X N L N q N 


(3.11) 


(3.12) 
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Here t , is the relative truncation error (or defect correction) 

N— 1 


between the grids so that the solution on the coarse grid is driven by 
the fine grid and the defect correction accounts for the difference in 
truncation error between the coarse and fine grids. 2 ® The analogous 
equation for (3.9) is given by 

L N-2 (Q N-2 ) = S N-2 + T N-2 (3.13) 

where 

(3.14) 


S = s 

N-2 N-1 N-1 


T = L (I N ~ 2 q C ) - I N " 2 (L (q C )) + I N " 2 T 
N-2 N-2 V N-1 4 N-1 } N-1 V N-1 VH N-1 ' ' N-1 N-1 


(3.15) 


Note that the relative truncation error on the N-2 grid is the sum of 
the relative truncation error between grids N and N-1 , as well as N-1 
and N-2. 


3.2 Algorithm for the Euler Equations 

For the steady Euler equations in generalized coordinates , equation 
(3.2) can be written as 

L (Q ) = 6“ + 6* F~ + 6“ G + + 6 + G - + 6~ H + + 6 + H~ = 0 (3.16) 

n' v n' 5 5 n n s c 

In the multigrid solution process, a forcing function arises on the 
coarse grids from restricting the residual equation on a fine mesh down 
to the coarser one. Since for the Euler equations, S = 0 in (3.2), the 
forcing function is the relative truncation error between the two 
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meshes. The resulting equation to be solved on any mesh can be 
written 

L.(Q.) = T i (3.17) 

where t. = 0 on the finest mesh and is the relative truncation error on 
l 

each of the coarser grids. The solution of equation (3.17) is updated 
by introducing a time derivative of the dependent variables to the left 
hand side so that the solution can be advanced in time using the 
approximate factorization methods described in Chapter Two. The 
resulting scheme written on mesh is given by 

M Aq? = -At(L.(q°) - T. ) = - AtR. (3.18) 

where M is the implicit operator of the scheme considered and LMq^) on 
the right hand side is due to the linearization of L^(Q^) from the 
backward Euler time integration. Note that even on the coarse meshes 
where x. is non-zero, equation (3.18) maintains the same form as the 
equation on the fine mesh. The result of this is that the coarse meshes 
can be updated using the same scheme as on the fine mesh with only a 
slight modification to the right hand side. 

There are several strategies for deciding when to switch from one 
grid level to another, generally falling under the categories of fixed 
or adaptive cycling algorithms. The strategy used in the present study 
is a fixed cycling strategy in which each global cycle consists of a set 
number of FAS cycles on each of the coarser grids. Recall that one FAS 
cycle consist of a smoothing step followed by coarse grid correction. 
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In addition, a predetermined number of iterations is performed on each 
grid level to smooth the errors before applying the coarse grid 
correction. 

The conserved variables are transferred to the next coarser grids 
each time by the rule 


Q. . = 1 Q. 

1-1 1 1 


( 3 . 19 ) 


where l| 1 is a volume weighted restriction operator which transfers 
values on the fine grid to the coarser one and is defined by 


= Ivq/Iv 


( 3 . 20 ) 


and the summations are taken over all the fine grid cells which make up 
the coarse grid cell. As shown in Appendix C, restriction of the 
dependent variables in this manner conserves mass, momentum, and energy 
in each of the cell volumes. The relative truncation error is 
calculated on the coarse grid as 


T = L (I?- ’q?) - 
1-1 1-1 11 11 


( 3 . 21 ) 


A i-1 

where 1^ is the restriction operator for the residual defined as 


if -1 R. = I R. 
11 ^ 1 


( 3 . 22 ) 


where again, the summation is over the cells on the fine grid which make 
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up the coarse grid. By summing the residuals, the surface integral of 
the fluxes crossing the cell boundaries on the coarse grid are the same 
as would arise by integrating around all the fine grid cells making up 
the coarse grid (see Appendix C). On this grid, several iterations of 
the approximate factorization scheme can be conducted to get an 
approximation to the steady solution on with the right hand side 

modified to include the relative truncation error. If only one coarse 
grid is used to correct the finest grid, the result is the simple FAS 
two-level cycle. On the other hand, if more grid levels are introduced 
so that one or more FAS cycles can be recursively carried out on each 
subsequent coarse grid level to get a better approximation to r then 

a multilevel algorithm results. When only one FAS cycle is carried out 
for each of the coarser grids, the resulting global cycling stragegy is 
termed a V-cycle and is depicted in figure 7. Another cycle of 
interest, which is shown in figure 8, is termed a W-cycle and results 
when two FAS cycles are used on each of the coarser meshes. Results 
will be shown in the next section using both types of cycles where on 
the coarsest mesh, three smoothing iterations are performed in lieu of 
solving the equations on the coarsest mesh exactly. The corrections on 
coarse meshes are passed to the next finer mesh using trilinear 
interpolation with no additional iteration steps between meshes. When a 
W-cycle is used however, note that an iteration is carried out at the 
beginning of each FAS cycle correction in order to smooth the high 
frequencies • 

In order to further clarify the multigrid procedure, the overall 
process is summarized below for an exemplary case where three grid 
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levels are used in a V-cycle. 

1 • Starting on the finest grid, smooth the errors by doing one 
iteration of equation (3*18) with 'T = 0. 

2. Calculate the residual on the fine grid from equation (3.3) where 
L N (q N c ) is given by equation (3.16). 

3. Restrict the dependent variables to the first coarse grid, , 

using equation (3.19). 

4. Restrict the residual from the finest grid to G n _.j using equation 
(3.22) and calculate the relative truncation error using (3.21). 

5. Calculate the right hand side of equation (3.18) and update the 
solution on mesh G N-1 * (This serves to smooth the errors on this 
grid so that a coarser grid can be introduced.) 

6. Calculate the residual on this mesh using equation (3.8). Note 
that this can be written as 

V, ■ L s-i (, S-.> - Vi < 3 - 23 ’ 

Since t has been previously calculated, the residual is easily 
N-1 

Q 

calculated by simply calculating L (q„ ) from the most current 

N-1 N-1 

values of the dependent variables on the mesh and 


subtracting t 


N-1 


7. Restrict the dependent variables on to G N _2 using equation 

(3.19). 


8. Restrict the residual from equation (3.23) to the N-2 grid and 

calculate x from equation (3.21). 

N-2 

9. Calculate the right hand side of equation (3.18) and update the 
solution on this mesh. Since this is the coarsest mesh used in the 
present example, three iterations of equation (3.18) are done to 
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get an approximation to Q N _ 2 * step, the right hand side is 

updated to use the most current values of the dependent variables 

in L _(qf, _)• Note that t 0 will not change. 

N— 2 N-2 N— 2 


10. Calculate the correction on this mesh 
c 


11 


V - q- - I N - 2 q C 
N-2 y N-2 N-1 h N-1 

Pass the correction to the next finest mesh using trilinear 

interpolation and update the solution 
c v c 


(3.24) 


q ' «- q v + I N-1 V 

q N-1 4 N-1 N-2 N-2 


12 


(3.25) 

Note that steps 5-11 make up one FAS cycle on grid N-1 where steps 

6-1 1 constitute a coarse grid correction. At this point, if a W- 

cycle were being employed, another FAS cycle (steps 5-11) would be 

c 

done to further update q 

N-1 

Calculate the correction on the N-1 mesh as 

„ c N-1 c 

V N-1 q N-1 “ Im q * 


N N 

13. Pass this correction to the finest mesh and update the solution 


c c N 

+ q« + 1 


v 

N-1 N-1 


N N 

14. Do one smoothing iteration using equation (3.18) to smooth the 


errors 
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Chapter IV 
RESULTS 

4.1 ONERA M6 

Three-dimensional subsonic and transonic flow computations over the 
ONERA M6 wing are shown below. Comparisons are made with experimental 
data at a Reynolds number of 11.7 million, 31 corresponding to conditions 
for which viscous effects are relatively small. The wing consists of 
symmetrical airfoil sections with a planform swept thirty degrees along 
the leading edge, an aspect ratio of 3.8, and a taper ratio of 0.56. 
Solutions are obtained for two mesh types, C-H and C-0, both of which 
are C-type mesh topologies around the airfoil profile. The C-H mesh, 
has uniform spacing in the spanwise direction whereas the C-0 mesh wraps 
around the wing tip, consequently leading to a more precise definition 
of the actual rounded tip geometry tested in the experiment. The C-0 
mesh has been generated with a trans-finite interpolationn procedure 
developed by Bruce Wedan of NASA Langley Research Center. The C-H mesh 
was obtained by simply stacking a series of two dimensional cross 
sections along the span. The surface mesh for both are shown in figure 
9 . 

The first computation is the ONERA M6 wing at transonic 
conditions: Mach number of 0.84 and an angle of attack of 3.06 

degrees. Figure 10 shows the effect of multigrid on the residual and 
lift history for a 193 x 33 x 33 C-H mesh, corresponding to 193 points 
along the airfoil and wake, 33 points approximately normal to the 
airfoil, and 33 points in the spanwise direction, 17 of which are on the 
wing planform. For this case, the Van Leer splittings are used with a 




(a) Residual history 


(b) Lift history 
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Figure 10. Effect of multigrid on convergence; ONERA M6; 

193 x 33 x 33 C-H mesh; M ro = 0.84; a = 3.06°. 
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V-cycle and four grid levels (a fine grid and three coarser ones). The 
multigrid scheme is very effective in accelerating convergence of both 
the residual and the lift. The residual is reduced to machine zero in 
400 cycles while the single grid scheme has only reduced the residual 
between one and two orders of magnitude. The benefit of multigrid is 
especially pronounced in the lift history where the final lift value is 
obtained to within 0.1 percent of its final value in 41 cycles. This is 
a dramatic improvement over the single grid result which required more 
than 400 iterations to settle in on a final lift coefficient. It should 
be noted that for all the cases considered, several cycles (usually 5) 
were run with first order spatial differencing before switching to 
second order. 

A comparison of convergence rates between the three schemes 
discussed in Chapter 2 is shown in figure 11 for identical conditions as 
given above with the exception that only every other point from the 193 
x 33 x 33 mesh is used, resulting in a 97 x 17 x 17 C-H mesh. For this 
size mesh, only two coarser grids are used. The 3-factor, spatially- 
split algorithm demonstrates a higher rate of convergence than either of 
the two factor schemes yielding a spectral radius of approximately 
0.898. The 2-factor scheme in which the implicit operator is split 
according to the sign of the eigenvalues displays the slowest 
convergence rate with a spectral radius of 0.93. It should be pointed 
out however that even though the spectral radius using this scheme is 
not as good as for the spatially split scheme, this still represents a 
good improvement over a corresponding single grid spectral radius of 
0.98. All the runs on the 97 x 17 x 17 meshes were made at a CFL 
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Cycles 


1 1 . Comparison of convergence rate for the 3 
schemes? ONERA M6 wing? 97 x 17 x 17 C-H 
mesh? M ot = 0.84? a = 3.06°. 
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number of 7. This was determined experimentally to be about optimum and 
agrees well with the CFL number for best smoothing predicted by the 
stability analysis* In 64 bit precision, the computational rate using a 
V-cycle and three grid levels for the 3-factor scheme is about 75 
microseconds per grid point per cycle whereas the 2-factor eigenvalue 
split and combination split schemes exhibit computational rates of 140 
and 85 microseconds per grid point per cycle, respectively. It should 
be noted that the computational rate is decreased by approximately 40 
percent when the computations are done in 32 bit precision. Due to the 
higher performance of the three-factor spatially split algorithm in both 
the convergence rate and the computational rate, it is used exclusively 
in the results that follow. 

The effect of using a W-cycle over the previously used V-cycle on 
the residual is shown in figure 12 for the 97x17x17 C-H mesh. An 
improvement using a W-cycle in the convergence rate is apparent. In 
addition, the lift coefficient is obtained to within 0.3 percent of its 
final value in only 14 cycles and to under 0.1 percent in 24 cycles. 

This is an improvement over the V-cycle which took 37 cycles to get the 
error in lift below 0.1 percent. Although the work involved for a W- 
cycle is more than for the V-cycle due to the extra smoothing iterations 
on the coarser grids, the time required per cycle only increased by 
about 13 percent over a V-cycle. Therefore, even though more work is 
involved for each cycle, a net gain is still achieved by employing the 
W-cycle. A summary of results for this case is given in table 2 for 193 
x 33 x 33 and 97 x 17 x 17 for both C-H and C-0 type grids. The table 
includes the spectral radius based on cycles and the number of cycles 






Mesh size and 
type of cycle 


97 x 17 x 17 C-H 
V-cycle 
W-cycle 


97 x 17 x 17 C-0 
V-cycle 
W-cycle 


193 x 33 x 33 C-H 

V-cycle 

W-cycle 


193 x 33 x 33 C-0 
V-cycle 
W-cycle 


Cycles required to obtain 


C£ to 0,3% C^ to 0.1% C^ to 5 
of final value of final value decimal places 


Spectral 

radius 


.898 

.871 


,912 

.879 


0.948 

0.923 


.952 

.926 


Table 2 Summary of results for ONERA M6 wing 
M„ = 0.84; a = 3.06°. 
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required to obtain the lift coefficient to within 0*3 and 0.1 percent of 
its final value as well as how many cycles were required to obtain the 
lift to five significant digits. Note that the number of cycles 
required for the W-cycle to obtain the lift coefficient is relatively 
insensitive to the number of grid points. 

Figure 13 shows the upper surface pressure distributions on the 193 
x 33 x 33 C-H mesh as well as the 193 x 33 x 33 C-0 mesh. The wing 
under these conditions exhibits both a swept shock emanating from the 
apex and a nearly normal shock emanating from the root, which coalesce 
at about 80 percent of the span to form a single shock. Figure 14 shows 
upper surface pressure contours. 

In figure 15, pressure coefficients obtained using the Van Leer 
splitting on both the 97 x 17 x 17 and 193 x 33 x 33 C-0 meshes are 

compared with experimental data at six spanwise locations. The 

computations are obtained at the same spanwise locations as the 
experimental data by linear interpolation. The computations on both 
meshes agree reasonably well with experiment for each spanwise location? 
the effect of the finer mesh is to resolve the leading edge suction 
pressures and shock positions, and improve the agreement with experiment 
at the most outboard station. Results obtained with the Steger-Warming 
splitting on the same two meshes are compared with experimental data in 
figure 16. The computations are nearly identical to the previous ones 
with small differences occurring near the shock regions. 

The next three-dimensional test case is the ONERA M6 wing at a 

freestream Mach number of 0.699 and an angle of attack of 3.06 

degrees. At these conditions, the flow remains subsonic over the entire 
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193 x 33 x 33 C-0 mesh 193 x 33 x 33 C-H mesh 


Figure 13. Upper surface variation of pressure coefficient; 
ONERA M6 wing; = 0.84; a = 3.06°. 
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wing. Results were obtained for this case on a 97 x 17 x 17 C-O mesh, a 
97 x 17 x 17 C-H mesh, and both a 193 x 33 x 33 C-H mesh and C-0 mesh. 
Figure 17 shows the residual history for both the 97 x 17 x 17 C-0 and 
C-H meshes using the Van Leer splittings and a V-cycle. The convergence 
rate on the C-H mesh is slightly better than on the C-0 mesh. Machine 

zero is reached for the C-H mesh in approximately 200 cycles while the 

C-0 mesh requires about 300 cycles, corresponding to an asymptotic 
spectral radius of .891 and .926, respectively. For both meshes, the 
lift was obtained to less than 0.1 percent of its final value in less 
than 28 cycles requiring only about 46 seconds of computer time. On the 
193 x 33 x 33 C-H mesh, a spectral radius of 0.929 was obtained with the 
multigrid algorithm while a spectral radius of 0.95 was obtained on the 
same sized C-0 mesh. When using a W-cycle, a spectral radius of 0.866 
is obtained for the 97 x 17 x 17 C-H mesh and one of 0.891 is obtained 

for the C-0 mesh. Using the 193 x 33 x 33 mesh, the spectral radius 

using the W-cycle is also about 0.89 for the C-H mesh and 0.912 for the 
C-0 mesh. A summary of results is given in table 3 similar to those 
shown in table 2. 

The pressure distributions on the 97 x 17 x 17 C-0 mesh and the 193 
x 33 x 33 C-H mesh are compared with experiment at six spanwise stations 
in figure 18 using the Van Leer splitting and in figure 19 using the 
Steger-Warming splittings. At the inboard stations, the results for 
both meshes are essentially identical and compare well with 
experiment. At the outboard station, however, the pressures computed on 
the C-0 mesh agree much closer to experiment due to the increased 


resolution at the tip 



Log 
(residual) 



Cycles 


97 x 17 x 17 C-H mesh 

97 x 17 x 17 C-0 mesh 


Figure 17. Multi grid convergence for C-H and C-0 meshes; 

ONERA M6 wing; M m = 0.699; a = 3.06°; V-cycle 


Cycles required to obtain 


Mesh size and 
type of cycle 

C£ to 0.3% 
of final value 

C£ to 0.1% 
of final value 

C£ to 5 

decimal places 

Spectral 

radius 

97 x 17 x 17 C-H 




0.891 

V-cycle 

19 

28 

38 

W-cycle 

11 

19 

33 

0.866 

97 x 17 x 17 C-0 





V-cycle 

21 

22 

48 

0.926 

W-cycle 

13 

15 

31 

0.891 

193 x 33 x 33 C-H 




HH 

V-cycle 

29 

37 

71 

Em 

W-cycle 

11 

21 

37 

»■ 

193 X 33 x 33 C-0 






V-cycle 

W-cycle 


0 

0 


21 

11 


38 

15 
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.950 

.912 
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4 . 2 Analytic Forebody 

The last test case considered is an analytically defined forebody 
for which experimental data is available at supersonic Mach numbers* 

The grid used, shown in figure 20 along with static density contours, 
was a 49 x 49 x 49 grid with a symmetry plane along the centerline. The 
conditions correspond to a freestream Mach number of 1.7 and an angle of 
attack of 0 degrees; this leads to an oblique shock at the nose and 
supersonic flow over the entire length of the body. 

The residual and lift history obtained using a V-cycle and the Van 
Leer splittings are shown in figure 21 . As can be seen, the residual is 
reduced 3 orders of magnitude in only 50 cycles (10 of which were first 
order accurate) and an asymptotic spectral radius of 0.83 is achieved 
based on the last thirty cycles. The lift is obtained to 0.3 percent 
error of the final value in only 22 cycles. The pressure distribution 
compared with experimental data at the forebody symmetry plane is shown 
in figure 22 over both the leeward and windward sides of the body. As 
seen, the pressure coefficients compare well with the experimental data 
over both the lower and upper surfaces of the body. 




(a) Residual history 



Cycles 


Figure 21 . Convergence history for ana 
forebody; M =1.7; a=0°. 


(b) Lift history 
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Figure 22. Calculated pressure coefficient comparison with 
experiment at forebody symmetry 
plane; = 1.7; a = 0°. 
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CHAPTER V 
CONCLUSIONS 

Multigrid acceleration has been applied to the three-dimensional 
flux-split Euler equations in generalized coordinates* Three implicit 
schemes have been used to smooth the errors at each grid level* Results 
from a linearized stability analysis of the coupled equations for each 
of the schemes agree well in overall trends with the numerical 
experiments and indicate that the 3-factor spatially-split algorithm is 
conditionally stable (up to a CFL of about 20) but offers a slightly 
better smoothing rate than the other two schemes and hence, the best 
multigrid performance* The stability analysis also indicates that the 
other two schemes, both two-factor schemes, are less sensitive to CFL 
variations. Results obtained for subsonic and transonic flow over the 
ONERA M6 wing and supersonic flow over an analytically defined forebody 
are compared with experimental data for cases with weak viscous 
effects. In the wing calculations, two methods of splitting the flux 
vectors were compared with experiment; the splitting of Steger-Warming 
and Van Leer. For both the subsonic and transonic cases, both methods 
of splitting the flux vector agree well with experiment and with each 
other. 

Results obtained for a series of subsonic, transonic, and 
supersonic flows demonstrate a substantial improvement in convergence 
rate using the multigrid algorithm in comparison to the single grid 
algorithm. Using a W-cycle, solutions can be obtained in as few as 19 
cycles for transonic conditions on a 193 x 33 x 33 mesh whereas a V- 
cycle takes about 41 cycles to reach the same level of accuracy (final 
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lift coefficient to within 0.1 percent). Using a V-cycle, a spectral 
radius of 0.891 and 0.898 is obtained for a 97 x 17 x 17 wing solution 
at subsonic and transonic conditions. A W-cycle for the same cases 
results in spectra radii of 0.866 for the subsonic case and 0.871 for 
the transonic case. In addition, the W-cycle is less sensitive to the 
size of the grid than the V-cycle for obtaining the final lift 
coefficient. Both of these save an order of magnitude in computing time 


over a single grid 
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APPENDIX A 

TRANSFORMATION TO GENERALIZED COORDINATES 
The three dimensional Euler equations in Cartesian coordinates and 
strong conservation law form are given by 


3Q 3F 3G 3h 

3t 3x 3y 3z 


( A. 1 ) 


where 


/ p 


pu 


f P v 1 


f ^ 

1 pu 


2 

pu + p 


puv 


paw 

< pv 


puv 

G = ' 

2 

pv + p 

H = ' 

pvw 

2 

I pw 

1 F = 

puw 


pvw 


pw + p 

' e 


(e+p)u ' 


1 (e+p)v ' 


1 (e+p)w ' 


12 2 2 

and p = (y_i)(e- — p(u + v + w )) 


Using the chain rule and the body-fitted coordinate system given by the 
steady transformation 


5 = 5(x,y,z ) n = n(x,y,z ) C = C(x,y,z) t = t 


(A. 3 ) 


the Euler equations can be recast as 


30 „ 3F 3F 3F , 3G 3G 3G 

3t + ? x 3C + n x 3n + ? x 3? + \ 85 + n y 8u + ? y 3? 


(A. 4) 


„ 3H 3H 3H 

4* E rv f. + T| “TT"" + £ fN _ — 

^z 3? z 3n z 35 


0 
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Now, again using the chain rule, the derivatives with respect 
to x, 5, r\, and ? can be written in matrix form as 


/ J_\ 

/ 1 o n o \ 

/4r\ 

3T 


1 


/ ^ 

9 




_3_ 

3? 


0 x 5 y ? 


3x 

_a 

= 



3 

3n i 


0 *» y n % 


i * 

, 3 / 


1 


3 i 

W 

\° X C y 5 Z C/ 


\ fe/ 


from which Cramer^ rule can be used to solve for the x,y,z and t 
derivatives* Using these to evaluate the metric terms gives 


- J(y n z ?- we 


5 y = j(z n x c - x n z c ) 


5 = J(x y - y x ) 

z n c n c 


-i 


where j = x f (y z - 


\ • J(z 5 y c" y f C 1 5 x 

^■•"VrVs 1 

V«> - Wc - Vc 1 


- J(y fn- z 5 y n> 

= JI Vn‘ x 5 z n ) 

" J<x ? y .r y S*n 1 

* 2 { (x n y c - Vt 1 


(A. 6) 


(A. 7) 


In order to regain the strong conservation law form, equation (A* 4) can 

now be multiplied through by J"" 1 and rearranged using the chain rule on 

certain terms* For example, the term j" 1 £ F can be rewritten as 

x ^ 


J "V« = -§c» J ' ,{ x> p)1 - 


(A. 8) 
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After rewriting all the appropriate terms and noting that many parts of 
the resulting equation can be shown to be zero by substituting equations 
(A. 6) for the metrics, the Euler equations can be written in generalized 
coordinates maintaining the strong conservation law form as 


9Q 9F 3G 8H 

TF + ‘55 +_ 5n +_ 5c" 0 

~ -1 

where Q = J Q 

* -1 

F = J (CF+SG+5H) 
x y z 

" -1 

G = J (riF+r|G + nH) 
x y z 

A -1 

H = J (5F+ ?GUH) 
X y z 


( A* 9 ) 


(A. 10) 


Using equations (A. 2) in (A. 10), the flux vectors can be further written 
in an alternate form as 


pu V 

pV 

PW 

puU + 5 X P J j 

pUV + T) P i 

X 1 J 

puW + ? x p 

PVU + ? p [ . \ 

l G = J 

pvv + j „ 1 | 

pvW + C^p 

pwu + 5 p 1 

Z 1 1 

pwv + n p 1 J 1 

1 z I 1 

pwW + ? z p 

(e + p)U ' ’ 

Me + p)V ' ' 

i (e + p)W 



where U, 

U = £ u 
x 

v = n u 

X 

W = c u 
x 
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V, and W are the contravariant velocities defined as 


+ £ v + 5 w 

y z 

+ n v + n w 

y z 

+ c v + c w (A. 12) 

y z 



71 


APPENDIX B 

SPLITTING THE FLUX VECTORS IN GENERALIZED COORDINATES 
The Van Leer method of splitting the flux vectors was originally 
given only for a Cartesian coordinate system.^ For example, the split 
flux vectors in the x-direction was given in terms of the local one- 
dimensional Mach number = u/a. For supersonic flow, i.e*, 

|mJ > 1 , we have 


F + = F, F" = 0, M > 1 


F + = 0 , F“ = F, M < -1 


(B.1) 


and for subsonic flow, m x < 1 


F ± 


±pa{l(M ± I)} 2 = f * 

2 x mass 


f * {( Y-1 )u ± 2a }/ y 

mass 1 ' 


f ± V 

mass 


f * w 
mass 


(B.2) 


f * c t{(Y-1)u ± 2a} 2 /{2(y 2 -1 ) } + (v 2 + w 2 )/2] 
mass 


For many applications, however, it is advantageous to construct 
generalized (body-fitted) coordinate systems of the type 
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5 = S(x,y,z) n = n(x,y,z) C = £(x,y,z) x = t (B.3) 

where, in the present work, the transformation is chosen so that the 
grid spacing in the computational domain is uniform and of unit 
length. In the discussion that follows, the superscript * indicates 
variables in generalized coordinates while superscript “ indicates 
variables in a locally Cartesian system. If no superscript is used, 
Cartesian coordinates are assumed. The strong conservation form of the 
Euler equations in generalized coordinates is given by 


A A 


9Q 3F 

TF + If 


A 

dG 


9h 


0 


( B. 4) 


For the purpose of determining a generalized splitting for F, only 
the derivatives in the 5 and t directions are considered while 
the ti and £ derivatives are treated as source terms. For determining 

A 

the splitting of F, equation (B.4) is transformed by a local rotation 

* 

matrix in order to decompose the flux vector F into components normal 
and tangential to a 5=constant cell face. The rotation matrix is given 
by 


T 


/’ 


0 

0 

0 

0 


x 


where 5,5, and 
x y 

a ^-constant line. 



!l 



(B.5) 


are components of a unit vector if normal to 

A A 

The t^ and r^ are components of vectors which are 
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normal to it* and to each other so that the three vectors form a locally 
Cartesian coordinate system. Note that an infinite number of vectors 
normal to ^ exists which form a locally Cartesian coordinate system. 
These vectors however are arbitrary and their exact specification is 
unnecessary. Multiplication of equation (B.4) with the matrix T then 
yields 


Q 


t 



-T 


G + T Q + F 
n t C 



(B. 6 ) 


where 



(B • 7 ) 


F 



grad £ 
J 


pu 

puu + p 

puv 

puw 

(e + p) u 


(B.8) 


The rotated velocity component u is the velocity normal to a line of 
constant £, representing the scaled contravariant velocity component, 
while v and w are normal to u and to each other 


u=(gu+5v+£w) 
x y z 


(B.9) 


v = (t u+tv + tw) 
x y z 


(B. 1 0 ) 
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AAA 

w = (ru+rv+rw) (B.11) 

x y z 


The transformed flux F is of the same functional form as the Cartesian 
flux vector and thus can be split according to any splitting developed 
for Cartesian coordinates. Therefore, equations for both the Steger- 
Warming and the Van Leer splittings can be used to split the flux 
vector F after replacing the Cartesian velocity components u, v and w by 
the rotated velocity components u, v and w. Applying the rotation T to 
equation (B.4) simply allows the flux vector to be split in a one- 
dimensional fashion, along a coordinate axis perpendicular to the cell 

A 

interface. After splitting F, the appropriate splitting for F is 
determined by applying the inverse transformation matrix T" 1 to equation 
(B.6), leading to 


Q +(F + F ) + G + H = 0 

t 5 n 5 


(B • 1 2 ) 


with 


p * = T' 1 p * . l gCadU) 


± 

mass 

± U (- u ± 2a) /Y + u] 

mass x 


f * o U, (- u ± 2a)/Y + V] 
mass y 

- ± 


f ± 

energy 


where 

f* = ± pa(M ± 1) 2 /4 
mass C 


(B.13) 
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= f maeet{-(Y - 1)U 2 ± 2(Y - 1)ua 
energy mass 


_ 2 . ,, 2 „ . ,2 2 2 . 

+ 2a }/(Y - 1) + (u + v + w )/2] 


Note that the inverse transformation restores the original form of the 
equations, i.e., no additional source terms arise and the form 

A A A A 

of G and H is unaffected. This allows a splitting of G and H similar to 

A 

the splitting of F shown above. 

Carried out with the Steger-Warming Cartesian splitting starting 
with equation (B.8) yields 


F = F + F + F 
1 2 3 

where 


- _ | grad ( g) [ «r Y-1 
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1 Y 1 pw 
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U 
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(B.16) 


(B.17) 


(B.18) 


Applying the inverse transformation to (B.14) gives 
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A A 


A ^ A 

F = T F = F 1 + f 2 + F 3 


(B.19) 
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pu 
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1 JY \ fw 
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p , 2 , 2 , 2. 

j (u + v + w ) 
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pu ± pa 5 

aA 

pv ± pa? 
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2,3 4,5 J2y I pw ± pa E 


(B.21) 


paU 

e + p ± -| — / 

grad ( 5) 


This is identical to the generalized splitting given in references 3 and 

7. 
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APPENDIX C 

RESTRICTION OPERATORS 


The dependent variables are transferred from a fine mesh to a 
coarse mesh so that the mass, momentum, and energy contained in the 
coarse grid cell is the same as that contained in the fine grid cells 
which compose the coarser grid. Referring to figure 23, for a two 
dimensional example, the mass, momentum and energy in any given cell is 
given by Q*V where V is the volume of the cell and Q is the column 
vector of dependent variable representing the conserved quantities per 
unit volume. Since cell A is comprised of the smaller cells a, b, c, 
and d, a relationship is easily established for transferring the 
dependent variables which conserves mass, momentum, and energy. 


q a - 


e a v a + 


Q b v b + 
V 


2c V c + 


2d V d 


(C.1) 


The restriction of the residuals is also guided by conservation 
laws. The steady Euler equations can be written in integral form as 



(C. 2) 


Here, the integral is the surface integral over the volume 

considered, F is the flux of mass, momentum, or energy across the 
* 

boundaries, and n is the outward pointing unit vector normal to the 
boundary. 

Considering the two dimensional case given above for simplicity, 
the integral around the large volume shown in figure 23 is given by 
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A 

-> 
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nds = (F . 

n) 1A S 1A + 

(F . 

n) 2A S 2A 

(C.3) 

, + 

A 

-► 



+ (F . 

n) 3A S 3A + 

(F . 

n) 4A S 4A 



Now, the integral along each of the larger faces is the sum of the two 
integrals along each of the smaller ones. For example 


(F 


n) 1A S 1A 


(E 1 . n) 1b S 1b + (F . n ) 1d s i d 


(C.4) 


Also, note that since the outward pointing normals on adjacent cell 
boundaries point in opposite directions, several terms which share a 
common boundary will cancel. For instance 


(F 


n) , S, 
1c 1c 


- (F • n) 3d S 3d 


(C.5) 


Therefore, by performing the integrations around each of the smaller 
cells and adding them together, it is seen that the integral around the 
larger cell is simply the sum of the integrals around each of the 
smaller ones. 


y A A 

/ F • nds = £ / F . nds (C.6) 

large cell small cells 

In order to better relate this specifically to the Euler equations, 
consider the steady continuity equation given by 

J pU • nds = 0 


(C. 7) 
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U = ui 


+ vj 


(C.8) 


n 


= k i + k 
x y 


j 



(C.9) 


A normal to a £=constant line is given by using k = 5 in (C.9) and the 
normal to an n=constant line is obtained using n in the same manner. 
Now, the length of each face is given by 



where k is again chosen to be 5 or r] depending on which face is desired 
and J represents the inverse of the cell volume (i.e. the Jacobian). 
Using equations (C.3) (with F = pU) , (C.7), (C.8), and (C.9), the 
integral around cell a in figure 23 can be written as 


/ ptf 


ndS = “ ( ZT>3a + 


(-PZ) - (PV) 

J '2a J ' 4a 


where 


(C.11) 


U = g u + 5 V 

x y 

v = n u + n v (c. i 2) 

x y 

When calculated numerically, this is simply the residual for the 
continuity equation. Similar results are obtained for the momentum and 
energy equations. Therefore, the residuals on a fine grid can be 
transferred to a coarser grid so that the integral of the fluxes over 
the cell boundary is conserved simply by summing the fine grid residuals 
which make up the coarse grid. 
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